from pyLFDA import LFDA
from IPython.display import SVG, display, Image
MDAnalysis : INFO MDAnalysis 2.0.0 STARTED logging to 'MDAnalysis.log' MDAnalysis : INFO MDAnalysis 2.0.0 STARTED logging to 'MDAnalysis.log' MDAnalysis : INFO MDAnalysis 2.0.0 STARTED logging to 'MDAnalysis.log'
the LFDA() class is central in running analysis on your simulations. \
\
Arguments:
experiment_name : Name of the experiment. Uses this to create a directory to store outputs in. If not specified time-stamp of experiment will be used.pdb_filename : Path of the PDB file to be used.gro_filename : Path of the GRO file to be used.trr_filename : Path of the TRR file to be used.tpr_filename : Path of the TPR file to be used.ndx_filename : Path of the NDX file to be used.gfda_version : Version of Gromacs FDA to be used. Creates a directory with the name to store it and uses it for further experiments. Currently, we support these cersions only - ['v2020.4-fda2.10.2', 'v2020.3-fda2.10.1', 'v2020.3-fda2.10', 'v2020-fda2.10', 'v2019.3-fda2.9.1', 'v2018.7-fda2.9.1']. #To do analysis, first we create an LFDA() object.
#It stores all files required to do the experiment.
#The argument names are self explanatory.
experiment = LFDA(experiment_name="lfda_het_cg", trr_filename="het-cg/step7_production.trr", tpr_filename="het-cg/step7_production.tpr", ndx_filename="het-cg/index.ndx", pdb_filename="het-cg/step7_production.pdb", gro_filename="het-cg/step7_production.gro", gfda_version="v2020.4-fda2.10.2")
#The specified version will automatically be downloaded and be used for any subsequent experiments if the files are run from the same directory.
Making MDA Universe from PDB and TRR file Parsing GRO file to calculate numbers of atoms, atoms information and box vectors /mnt/c/Users/bhava/Desktop/work/lfda/pyLFDA/het_memb/het-cg/step7_production.gro file parsed. with 14564 atoms in 0.07831740379333496 seconds
run_fda() create .pfi file and then generating a PFA file using GROMACS FDA. \
\
Arguments :
group1 : 1st group selectedgroup2 : 2nd group selectedresidue_list : [group1, group2]pfi_filename : Name of the PFI file to be generated. It is inferred from the experiment class if None.pfa_filename : Name of the PFA file to be generated. It is inferred from the experiment class if None.#Run GROMACS FDA to calculate pairwise forces.
#This command creates a .pfa file containing pairwise forces between atoms.
experiment.run_fda(group1="POPS", group2="POPC", residue_list=["POPS", "POPC"], pfa_filename="pfa_file.pfa")
Creating PFI file /mnt/c/Users/bhava/Desktop/work/lfda/pyLFDA/het_memb/lfda_het_cg/pfi_Dec-17-2021_1935.pfi file created in 0.004339694976806641 seconds Running Gromacs FDA ResidueRenumber: auto Vector2Scalar: norm Pairwise interactions selected: all Pairwise forces for groups: POPS and POPC Binary mode: 0 Threshold: 1e-10 Normalize punctual stress per residue: 0 Ignore missing potentials: 1
:-) GROMACS - gmx mdrun, 2020.4-dev-20201029-4b5790511-unknown (-:
GROMACS is written by:
Emile Apol Rossen Apostolov Paul Bauer Herman J.C. Berendsen
Par Bjelkmar Christian Blau Viacheslav Bolnykh Kevin Boyd
Aldert van Buuren Rudi van Drunen Anton Feenstra Alan Gray
Gerrit Groenhof Anca Hamuraru Vincent Hindriksen M. Eric Irrgang
Aleksei Iupinov Christoph Junghans Joe Jordan Dimitrios Karkoulis
Peter Kasson Jiri Kraus Carsten Kutzner Per Larsson
Justin A. Lemkul Viveca Lindahl Magnus Lundborg Erik Marklund
Pascal Merz Pieter Meulenhoff Teemu Murtola Szilard Pall
Sander Pronk Roland Schulz Michael Shirts Alexey Shvetsov
Alfons Sijbers Peter Tieleman Jon Vincent Teemu Virolainen
Christian Wennberg Maarten Wolf Artem Zhmurov
and the project leaders:
Mark Abraham, Berk Hess, Erik Lindahl, and David van der Spoel
Copyright (c) 1991-2000, University of Groningen, The Netherlands.
Copyright (c) 2001-2019, The GROMACS development team at
Uppsala University, Stockholm University and
the Royal Institute of Technology, Sweden.
check out http://www.gromacs.org for more information.
GROMACS is free software; you can redistribute it and/or modify it
under the terms of the GNU Lesser General Public License
as published by the Free Software Foundation; either version 2.1
of the License, or (at your option) any later version.
GROMACS: gmx mdrun, version 2020.4-dev-20201029-4b5790511-unknown
Executable: /mnt/c/Users/bhava/Desktop/work/lfda/pyLFDA/het_memb/v2020.4-fda2.10.2/bin/./gmx_fda
Data prefix: /mnt/c/Users/bhava/Desktop/work/lfda/pyLFDA/het_memb/v2020.4-fda2.10.2
Working dir: /mnt/c/Users/bhava/Desktop/work/lfda/pyLFDA/het_memb
Command line:
gmx_fda mdrun -rerun /mnt/c/Users/bhava/Desktop/work/lfda/pyLFDA/het_memb/het-cg/step7_production.trr -s /mnt/c/Users/bhava/Desktop/work/lfda/pyLFDA/het_memb/het-cg/step7_production.tpr -pfi /mnt/c/Users/bhava/Desktop/work/lfda/pyLFDA/het_memb/lfda_het_cg/pfi_Dec-17-2021_1935.pfi -nt 1 -pfa /mnt/c/Users/bhava/Desktop/work/lfda/pyLFDA/het_memb/lfda_het_cg/pfa_file.pfa -pfn /mnt/c/Users/bhava/Desktop/work/lfda/pyLFDA/het_memb/het-cg/index.ndx
Back Off! I just backed up md.log to ./#md.log.3#
Compiled SIMD: None, but for this host/run AVX2_256 might be better (see log).
The current CPU can measure timings more accurately than the code in
gmx mdrun was configured to use. This might affect your simulation
speed as accurate timings are needed for load-balancing.
Please consider rebuilding gmx mdrun with the GMX_USE_RDTSCP=ON CMake option.
Reading file /mnt/c/Users/bhava/Desktop/work/lfda/pyLFDA/het_memb/het-cg/step7_production.tpr, VERSION 2020 (single precision)
Residue number collision detected, residues will be renumbered.
Changing nstlist from 20 to 25, rlist from 1.212 to 1.26
Using 1 MPI thread
Non-default thread affinity set, disabling internal thread affinity
Using 1 OpenMP thread
WARNING: Using the slow plain C kernels. This should
not happen during routine usage on supported platforms.
Back Off! I just backed up traj.trr to ./#traj.trr.3#
Back Off! I just backed up ener.edr to ./#ener.edr.3#
WARNING: This run will generate roughly 5006 Mb of data
starting md rerun 'Martini system', reading coordinates from input trajectory '/mnt/c/Users/bhava/Desktop/work/lfda/pyLFDA/het_memb/het-cg/step7_production.trr'
trr version: GMX_trn_file (single precision)
Reading frame 9000 time 900000.000
GMX RUN completed in 615.8366215229034 seconds /mnt/c/Users/bhava/Desktop/work/lfda/pyLFDA/het_memb/lfda_het_cg/pfa_file.pfa PFA file is generated in /mnt/c/Users/bhava/Desktop/work/lfda/pyLFDA/het_memb/lfda_het_cg
Last frame 10000 time 1000000.000
NOTE: 15 % of the run time was spent in pair search,
you might want to increase nstlist (this has no effect on accuracy)
Core t (s) Wall t (s) (%)
Time: 615.627 615.642 100.0
(ns/day) (hour/ns)
Performance: 140341.243 0.000
GROMACS reminds you: "What is a Unix or Linux sysadmin's favourite hangout place? Foo Bar." (Anonymous)
#Running GROMACS FDA can be a time consuming process.
#If you want to continue experimenting on a .pfa file you had previously generated.
#you can just load it and continue working.
#you dont need to run `run_fda()` if you load your .pfa file.
experiment.load_pfa(pfa_filename="lfda_het_cg/pfa_file.pfa", group1="POPS", group2="POPC", residue_list=["POPS", "POPC"])
Loading PFA file generated by Gromacs FDA
The .pfa file generated by GROMACS FDA is parsed for further analyses. The file can be parsed in two ways depending on your use case.
Average Parsing - Allows for calculation of averaged properties over all frames. set mode argument to average.Framewise Parsing - Allows for calculation of either properties for a specific frame or for a moving window over all frames. set mode argument to framewise.#The .pfa file created, is parsed as average or framewise.
experiment.framewise=False #we specify that we need to parse pfa as average.
experiment.parse_pfa(file_name="generated_pfa_average.pfa")
Parsing PFA file as average Parsed PFA file /mnt/c/Users/bhava/Desktop/work/lfda/pyLFDA/het_memb/lfda_het_cg/generated_pfa_average.pfa created with 14564 atoms in 263.49041056632996 seconds Summed PFA file parsed in 0.06467533111572266 seconds
#Parsing .pfa can also be time consuming.
#You can just load the parsed .pfa to continue experimenting.
#The generated pfa file should be parsed as average to further parse it as average and similarly for framewise.
#using this function in "average" mode automatically sets "experiment.framewise" to "False".
experiment.parse_parsed_pfa(file_name="lfda_het_cg/generated_pfa_average.pfa", mode="average", group1="POPS", group2="POPC", residue_list=["POPS","POPC"])
Loading Average parsed PFA file Summed PFA file parsed in 0.029172658920288086 seconds
This function generates curvature plots using MDAnalysis and displays them along with pairwise forces of the selected groups and the angle they make with the z-axis. \ \ Arguments:
specific_frame: (int) Frame to calculate forces for. default: None.window: (int) Moving Window size. default: None.num_x_bins: (int) Number of bins in x-direction. default: 10.num_y_bins: (int) Number of bins in y-direction. default: 10.split: (bool) Split Calculations into Upper and Lower Membranes. default: False#Generate curavture plots
experiment.curvature()
MDAnalysis.analysis.base: INFO Choosing frames to analyze MDAnalysis.analysis.base: INFO Choosing frames to analyze MDAnalysis.analysis.base: INFO Starting preparation MDAnalysis.analysis.base: INFO Starting preparation
Creating curvature plot
MDAnalysis.analysis.base: INFO Finishing up MDAnalysis.analysis.base: INFO Finishing up
#visualise curvature plots
display(SVG("lfda_het_cg/curvature_averaged_10_10_P.svg"))
average_z_surface, average_mean_curvature and average_gaussian_curvature respectively.#Splitting the curavture into upper and lower membrane.
experiment.curvature(split=True)
[68.64947455 67.87826523 43.27094985] 1
MDAnalysis.analysis.base: INFO Choosing frames to analyze MDAnalysis.analysis.base: INFO Choosing frames to analyze MDAnalysis.analysis.base: INFO Starting preparation MDAnalysis.analysis.base: INFO Starting preparation
Creating curvature plot
MDAnalysis.analysis.base: INFO Finishing up MDAnalysis.analysis.base: INFO Finishing up MDAnalysis.analysis.base: INFO Choosing frames to analyze MDAnalysis.analysis.base: INFO Choosing frames to analyze MDAnalysis.analysis.base: INFO Starting preparation MDAnalysis.analysis.base: INFO Starting preparation MDAnalysis.analysis.base: INFO Finishing up MDAnalysis.analysis.base: INFO Finishing up
display(SVG("lfda_het_cg/curvature_averaged_10_10_P_Upper.svg"))
display(SVG("lfda_het_cg/curvature_averaged_10_10_P_Lower.svg"))
This function clusters the selected residues. \ \ Arguments:
attached_ligands: (str) Ligand 1 to be clustered. default: group 1 previously selected("POPS" in this example).lipids_to_cluster: (str) Ligand 2 to be clustered. default: group 2 previously selected("POPC" in this example).protein_residue_names: (list) All residues except for attached_ligands and lipids_to_cluster.box_side_length: (int) Length of a box size. default: 6.experiment.cluster(lipids_to_cluster="POPS", attached_ligands="POPC")
Making clutering plots
display(SVG("lfda_het_cg/cluster_POPS_POPC.svg"))
experiment.cluster(lipids_to_cluster="POPS", mode="single")
Making clutering plots
display(SVG("lfda_het_cg/cluster_POPS.svg"))
This function calculates the pairwise force between selected groups 1 and 2. \ \ Arguments:
specific_frame: (int) Frame to calculate forces for. default: None.window: (int) Moving Window size. default: None.experiment.force_graph()
Creating average force plot Average Force plots created and saved
display(SVG("lfda_het_cg/force_averaged_['POPS', 'POPC'].svg"))
This function calculates MSD values and calculates the diffusion coefficient using MDAnalyses.
experiment.msd()
MDAnalysis.analysis.base: INFO Choosing frames to analyze MDAnalysis.analysis.base: INFO Choosing frames to analyze MDAnalysis.analysis.base: INFO Starting preparation MDAnalysis.analysis.base: INFO Starting preparation
Calculating diffusion coefficient
MDAnalysis.analysis.base: INFO Finishing up MDAnalysis.analysis.base: INFO Finishing up
Diffution coefficient and MSD plotted
display(SVG("lfda_het_cg/MSD.svg"))
This section illustrates how to generate framewise and moving window plots.
#The .pfa file created, is parsed asframewise.
experiment.framewise = True
experiment.parse_pfa(file_name="generated_pfa_framewise.pfa")
Parsing PFA file as framewise Parsed PFA file /mnt/c/Users/bhava/Desktop/work/lfda/pyLFDA/het_memb/lfda_het_cg/generated_pfa_framewise.pfa created with 14564 atoms in 760.5916864871979 seconds PFA file parsed in 654.4875428676605 seconds
#using this function in "framewise" mode automatically sets "experiment.framewise" to "True".
experiment.parse_parsed_pfa(file_name="lfda_het_cg/generated_pfa_framewise.pfa", mode="framewise", group1="POPS", group2="POPC", residue_list=["POPS","POPC"])
Loading Framewise parsed PFA file Summed PFA file parsed in 965.686830997467 seconds
#create curvature plots for the 50th frame
experiment.curvature(specific_frame=50)
Creating 50 specific frame curvature plot
MDAnalysis.analysis.base: INFO Choosing frames to analyze MDAnalysis.analysis.base: INFO Choosing frames to analyze INFO:MDAnalysis.analysis.base:Choosing frames to analyze MDAnalysis.analysis.base: INFO Starting preparation MDAnalysis.analysis.base: INFO Starting preparation INFO:MDAnalysis.analysis.base:Starting preparation MDAnalysis.analysis.base: INFO Finishing up MDAnalysis.analysis.base: INFO Finishing up INFO:MDAnalysis.analysis.base:Finishing up
#visualise the plots
display(SVG("lfda_het_cg/curvature_framewise_10_10_P_50.svg"))
#Framewise lipid angles
display(SVG("lfda_het_cg/curvature_framewise_angles_framewise.svg"))
#create curvature plots for a moving window of size 20. Plots created are for the windows 0-1000, 1000-2000 and so on till 10k for our case.
experiment.curvature(window=1000)
Creating 1000 window size curvature plots
MDAnalysis.analysis.base: INFO Choosing frames to analyze MDAnalysis.analysis.base: INFO Choosing frames to analyze MDAnalysis.analysis.base: INFO Starting preparation MDAnalysis.analysis.base: INFO Starting preparation MDAnalysis.analysis.base: INFO Finishing up MDAnalysis.analysis.base: INFO Finishing up
display(SVG("lfda_het_cg/curvature_moving_window_10_10_P_window_5000_6000.svg"))
#Gangle isnt working
display(SVG("lfda_het_cg/curvature_moving_window_angles_windowed.svg"))
#create force plots for the 50th frame
experiment.force_graph(specific_frame=50)
Creating framewise average force plot Force plots created and saved for frame - 50
#visualise the plots
display(SVG("lfda_het_cg/force_specific_frame_POPS_50.svg"))
display(SVG("lfda_het_cg/force_specific_frame_POPC_50.svg"))
#create force plots for a moving window of size 20.
experiment.force_graph(window=20)
Creating framewise average force plot Force plots created and saved for moving window of size - 20
display(SVG("lfda_het_cg/force_moving_window_POPC_40_to_60.svg"))
display(SVG("lfda_het_cg/force_moving_window_POPS_40_to_60.svg"))
#Create new PDB files with bFactor Loaded.
experiment.bfactor_pdb(bfactor_pdb_filename="bfactor_pdb_combined", mode="combined")
Loading a new PDB file with bFactor PDB with BFactor values created in 1466.0400023460388 seconds
experiment.bfactor_pdb(bfactor_pdb_filename="bfactor_pdb_atomwise", mode="atomwise")
Loading a new PDB file with bFactor PDB with BFactor values created in 2111.4183835983276 seconds
This function calculates the angle between the vector of the selected lipids and the z-axis. The vector is defined between the P atom and the specified atom. The function allows you to calculate the angle individually for each lipid or group them together to calculate the average angle. \ \ Arguments:
selection: (list) List of lipids for which the angle has to be calculated.grouping: ('combine'/'individual') Whether to group the selected lipids or keep them individual.c_atom_name: (str) Name of the atom to which the lipid vector is to be defined. default: 'C4B'.experiment.angles(selection=['POPC', 'POPS', 'POPE', 'POSM'], grouping = 'individual', c_atom_name = 'C4B')
Calculating lipid angles with vector as P -> C4B
display(SVG("lfda_het_cg/angles_framewise_P_to_C4B_POPC_POPS_POPE_POSM.svg"))
experiment.angles(selection=['POPC', 'POPS', 'POPE', 'POSM'], grouping = 'individual', c_atom_name = 'C4B', split = True)
Calculating lipid angles with vector as P -> C4B
display(SVG("lfda_het_cg/angles_framewise_P_to_C4B_POPC_POPS_POPE_POSM_Upper.svg"))
display(SVG("lfda_het_cg/angles_framewise_P_to_C4B_POPC_POPS_POPE_POSM_Lower.svg"))
experiment.angles(selection=['POPC', 'POPS', 'POPE', 'PSM'], grouping = 'combine', c_atom_name = 'C4B')
Calculating lipid angles with vector as P -> C4B
display(SVG("lfda_het_cg/angles_framewise_combined_P_to_C4B.svg"))